/**
 * @brief    平面九节点四边形截面属性计算单元测试
 * 
 * @version  0.1
 * @author   Liu Xing
 * @date     2025.10.17
 * @copyright Copyright (c) 2025
 */

#include "DataStructure/Common/Constant.h"
#include "gtest/gtest.h"
#include "../Quad9Element.h"

using SG::DataStructure::Common::TOLERANCE;
using SG::Algebra::Matrixd;
using SG::FEMSolver::BeamElement::Quad9Element;
using SG::FEMSolver::BeamElement::EleGeomProp;

static Matrixd coord(9, 2, {
                            1.0, 2.0,
                            7.0, 2.0,
                            7.0, 6.0,
                            1.0, 6.0,
                            4.0, 2.0,
                            7.0, 4.0,
                            4.0, 6.0,
                            1.0, 4.0,
                            4.0, 4.0});


// 测试面积计算
TEST(Quad9ElementTest, CalculateGeomPropTest)
{
    Quad9Element ele;
    ele.Initialize();
    ele.CalculateGeomProp(coord);
    EleGeomProp eleGeoProp = ele.GetGeomProp();
    EXPECT_NEAR(eleGeoProp.Area, 24.0, TOLERANCE);
    EXPECT_NEAR(eleGeoProp.Qy, 96.0, TOLERANCE);
    EXPECT_NEAR(eleGeoProp.Qz, 96.0, TOLERANCE);
    EXPECT_NEAR(eleGeoProp.Iy, 416.0, TOLERANCE);
    EXPECT_NEAR(eleGeoProp.Iz, 456.0, TOLERANCE);
    EXPECT_NEAR(eleGeoProp.Iyz, 384.0, TOLERANCE);
}

// 测试计算单元积分项ComputeShearCoeffIntegralTerm()
TEST(Quad9ElementTest, ComputeShearCoeffIntegralTermTest)
{
    Quad9Element ele;
    ele.Initialize();

    // 构造线性分布 psi, phi
    Matrixd psi(9, 1), phi(9, 1);
    for (size_t i = 0; i < 9; ++i)
    {
        psi(i, 0) = static_cast<Real>(i);
        phi(i, 0) = static_cast<Real>(i) * 0.5;
    }

    const auto result = ele.ComputeShearCoeff(coord, {1.0, 1.0, 0.0}, 0.0, psi, phi);
    ASSERT_TRUE (result == Matrixd(3, 1, {133.34444444444429, 33.336111111111073,66.672222222222146}));

}

// 测试计算单元矩阵ComputeEleMatrix()
TEST(Quad9ElementTest, ComputeEleMatrixTest)
{
    Quad9Element ele;
    ele.Initialize();

    Matrixd Ke;
    std::vector<Matrixd> F;
    ele.ComputeEleMatrix(coord, {1.0, 1.0, 0.0}, 0.0, Ke, F);
    
    // Check
    ASSERT_TRUE (Ke(0, 0) == 0.6740740740740726);
    ASSERT_TRUE (Ke(2, 3) == -0.087037037037037204);
    ASSERT_TRUE (Ke(8, 8) == 6.1629629629629603);
    ASSERT_TRUE (F[0](6,0) == -16.000000000000007);
    ASSERT_TRUE (F[1](3,0) == 1.3333333333333295);
    ASSERT_TRUE (F[2](7,0) == 21.333333333333314);
}

// 测试计算单元扇形积
TEST (Quad9ElementTest, ComputeSectorProductGaussIntegral)
{
    Matrixd eleNodeOmega (9, 1, {2.0, 4.0, 6.0, 8.0, 16.0, 24.0, 14.0, 28.0, 42.0});

    // Test data
    Quad9Element quad9;
    quad9.Initialize ();
    auto IyzwTest = quad9.ComputeSectorProduct (coord, eleNodeOmega);

    // Target data
    Matrixd IyzwTarget (1, 2, {2688.0, 2720.0});

    // Check
    ASSERT_TRUE (IyzwTarget == IyzwTest);
}